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Abstract. Symplectic integration algorithms have become popular in 
recent years in long-term orbital integrations because these algorithms en- 
force certain conservation laws that are intrinsic to Hamiltonian systems. 
For problems with large variations in timescale, it is desirable to use a 
variable timestep. However, naively varying the timestep destroys the 
desirable properties of symplectic integrators. We discuss briefly the idea 
that choosing the timestep in a time symmetric manner can improve the 
performance of variable timestep integrators. Then we present a symplec- 
tic integrator which is based on decomposing the force into components 
and applying the component forces with different timesteps. This multi- 
ple timescale symplectic integrator has all the desirable properties of the 
constant timestep symplectic integrators. 



1. Symplectic Integrators 

Long-term numerical integrations play an important role in our understanding of 
the dynamical evolution of many astrophysical systems (see, e.g., Duncan, these 
proceedings, for a review of solar-system integrations). An essential tool for 
long-term integrations is a fast and accurate integration algorithm. Symplectic 
integration algorithms (SIAs) have become popular in recent years because the 
Newtonian gravitational A-body problem is a Hamiltonian problem and SIAs 
enforce certain conservation laws that are intrinsic to Hamiltonian systems (see 
Sanz-Serna & Calvo 1994 for a general introduction to SIAs). 

For an autonomous Hamiltonian system, the equations of motion are 

dw/dt = {w,H}, (1) 

where H(w) is the explicitly time-independent Hamiltonian, w = (q,p) are the 
2d canonical phase-space coordinates, { , } is the Poisson bracket, and d(= 3N) 
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is the number of degrees of freedom. The formal solution of Eq.(|l]) is 

w(t) = exp(i{ ,H})w(0). (2) 



If the Hamiltonian H has the form Ha + Hg, where Ha and Hb are sepa- 
rately integrable, we can devise a SIA of constant timestep r by approximating 
exp(r{ , H}) as a composition of terms like exp(r{ ,Ha}) and exp(r{ ,Hb})- 
For example, a second-order SIA is 

exp(^{ , H A }) exp(r{ , Hb}) exp(^{ , H A }). (3) 

For the gravitational iV-body problem, we can write H = T(p) + V(q), 
where T(p) and V(q) are the kinetic and potential energies. Then the second- 
order SIA Eq.(|3|) becomes 

P n+ i = Pn + 7^(9n)> Qn+l =Qn + Pn+1 = P„+I + ^(Qn+l)' 

(4) 

where F = —dV/dq and v = dT/dp; this is the familiar leapfrog integrator. For 
solar-system type integrations, a central body (the Sun) is much more massive 
than the other bodies in the system, and it is better to write H = H^ep + -f^int, 
where i^Kep is the part of the Hamiltonian that describes the Keplerian motion 
of the bodies around the central body and H{ nt is the part that describes the 
perturbation of the bodies on one another. Symplectic integrators using this 
decomposition of the Hamiltonian were introduced by Wisdom & Holman (1991), 
and they are commonly called mixed variable symplectic (MVS) integrators. 
The constant timestep SIAs have the following desirable properties: 

(1) As their names imply, SIAs are symplectic, i.e., they preserve dp A dq. 

(2) For sufficiently small r, SIAs solve almost exactly a nearby "surrogate" 
autonomous Hamiltonian problem with H = H + H eiI . For example, the second- 
order SIA in Eq.(|3|) has 

#crr = ^{{Ha, H b }, H b + \h a } + 0(r 4 ). (5) 

Consequently, we expect that the energy error is bounded and the position (or 
phase) error grows linearly. 

(3) Many SIAs (e.g., Eq.||) are time reversible. Note, however, that there are 
algorithms that are symplectic but not reversible. 

Fig. 1 shows the energy error AE/E of an integration of the e = 0.5 Kepler 
orbit using the constant timestep leapfrog integrator Eq.(||). (In this and all 
subsequent Figures, the orbits are initially at the apocenter r a = a(l + e), where 
a and e are the semi-major axis and eccentricity.) Fig. 1 illustrates that there 
is no secular drift in AE/E. 

2. Simple Variable Timestep 

For problems with large variations in timescale (due to close encounters or high 
eccentricities), it is desirable to use a variable timestep. A common practice 
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Figure 1. Energy error AE/E of an integration of the Kepler prob- 
lem using the constant timestep leapfrog integrator Eq.(Q). The eccen- 
tricity e = 0.5, the timestep r = P/4000, and P is the orbital period. 
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Figure 2. (a) Same as Fig. 1, but using the leapfrog integrator Eq.(Q) 
with the simple variable timestep r = To(r/r a ) 3 / 2 , where tq = P/2000. 
The errors for two neighboring initial conditions are also shown (see 
text), (b) Same as (a), but using the symmetrized timestep criterion. 




is to set the timestep using the phase-space coordinates w n at the beginning 
of the timestep: r = h(w n ). If the SIAs discussed in § 1 are implemented 
with this simple variable timestep scheme, they are still symplectic if we assume 
that the sequence of timesteps determined for a particular initial condition are 
also used to integrate neighboring initial conditions (see Skeel & Gear 1992 for 
another point of view). However, tests have shown that this and similar simple 
variable timestep schemes [and also r = h(t)] destroy the desirable properties 
of the integrators (e.g., Gladman, Duncan, & Candy 1991; Calvo & Sanz-Serna 
1993). In Fig. 2a we show the energy error of an integration with e = 0.5 and 
h(r) = To(r/r a ) 3 / 2 . Although the error is initially smaller than that in Fig. 1 (the 
integrations shown in Figs. 1 and 2 use nearly the same number of timesteps 
per orbit), it shows a linear drift. In Fig. 2a we also show the errors for two 
neighboring initial conditions (e = 0.5 ± 0.005). They were integrated using the 
timesteps determined for the e = 0.5 orbit. Note that these errors grow even 
faster. 



The degradation in performance is due to the fact that the properties (2) 
and (3) listed in § 1 are no longer true. The algorithm is not time reversible 
because the timestep depends only on the coordinates at the beginning of the 
timestep. Since H err depends explicitly on the timestep r (see, e.g., Eq.Jg]), a 
variable timestep changes the surrogate Hamiltonian problem that the integrator 
is solving from step to step. Thus the solution from t = to t n after n steps is 
not in general the solution of a nearby autonomous Hamiltonian problem. 

3. Time Symmetrization 

Hut, Makino, & McMillan (1995; see also Funato et al. 1996; Hut, these pro- 
ceedings) pointed out that the performance of a variable timestep integrator can 
be improved if time reversibility is restored by choosing the timestep in a time 
symmetric manner: r = r(w n , w n +\) with r(w n , w n +i) = r(w n +i,w n ). For 
example, r = [h(w n ) + h{w Tl+ \)\/2. We have tested this idea in detail by (i) 
integrating a series of problems (the pendulum, Kepler orbits, and the restricted 
3-body problem) using a symmetrized variable timestep leapfrog integrator and 
(ii) integrating the restricted 3-body problem using a symmetrized MVS inte- 
grator. We have found that time symmetrization usually reduces the drift in 
energy (or the Jacobi constant) to negligible levels. Fig. 2b shows an integration 
of the e = 0.5 Kepler orbit using the symmetrized leapfrog integrator, and it 
should be compared to Fig. 2 a. 

In Fig. 2b we also show the errors for two neighboring initial conditions that 
were integrated using the timesteps determined for the e = 0.5 orbit. There 
is no improvement in the error growth for the neighboring initial conditions. 
Integrating an orbit using the timesteps determined for another orbit may seem 
to be somewhat artificial, but similar situations do occur in realistic integrations. 
For example, if we integrate an iV-body system using a shared timestep scheme, 
the timestep used by all of the particles is determined by the few having the 
strongest close encounter. In practice, since the symmetrized timestep criterion 
depends on io n +i, a symmetrized integrator also has the disadvantage that it 
requires iteration and can be significantly slower than the original integrator 
(unless the original integrator is also implicit). 

4. Multiple Timescale Symplectic Integrators 

In this section we describe a "variable timestep" SIA that has all the desirable 
properties of the constant timestep SIAs. The algorithm is based on ideas pro- 
posed by Skeel & Biesiadecki (1994; see also MacEvoy & Scovel 1994). It is also 
similar to the individual timestep scheme of Saha & Tremaine (1994). 

For simplicity, let us consider in particular the Kepler problem with T(p) = 
\p\ 2 /2 and V(r) = —l/r. We choose a set of cutoff radii n > r2 > • • • and de- 
compose the potential V into Vi, or equivalently the force F into Fi = —dVi/dr, 
such that (i) F = J2iZoFi, (ii) F{ (except Fq) is zero at r > rj, and (iii) F{ 
is "softer" than F; l+ \. The force Fi is to be applied with a timestep r». If we 
assume that Tj/rj+i = Mj + i is an integer, we can apply the second-order SIA in 
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Figure 3. Decomposition of the force F(r) = 1/r 2 into F c i or F p ^. 
The i = 0, 1, and 2 components, with r\ = 1 and rj/rj + i = \/2, are 
shown. 



Eq.(||) recursively to obtain the following second-order algorithm: 

exp(r { , H}) « exp(^Ko) exp[r (Z) + #1 + tf 2 + • • •)] exp(^K ) 



exp(yKo) 



exp(^tfi) exp[n(£) + K 2 + K 3 + • • •)] 



x exp(^-iTi) 



Afi 



(6) 



where D = { , T} and K j = { , V^}. Hereafter, we adopt Mj = 2. 

The multiple timescale algorithm Eq.(|6|) has an overall timestep To, but it is 
effectively a variable timestep scheme because the recursion terminates at level 
i if Ki + \ + Ki + 2 + • • • = during a substep of length t%. For example, if the 
particle is in the region r% > r > r 2 during an overall step, Eq.(||) reduces to 

eM^K + T ^K X ) exp(^D) exp(^tfi) exp(^D) exp(^ET 1 + ^K ). (7) 

Since Eq. (|6|) is based on the recursive application of Eq. (|3|) , it has all the desir- 
able properties of a constant timestep SIA. It is obviously symplectic and time 
reversible. We can also derive the surrogate autonomous Hamiltonian solved by 
this integrator. For example, if we decompose V into two levels Vq and V± only, 
the error Hamiltonian is (cf. Eq.||) 

+ O(r 4 ). (8) 

After some experiments, we found two force decompositions that work well. 
If we write F c ^ = F c ,o and F C) i = F c ^ — F c ^-\ for i ^ 0, one of the decompo- 
sitions is 



c,i—l 



if r > rj, 

9 (r/rj) 2 — 5 (r/rj) 6 r/Arf if r < rj. 



(9) 
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Figure 4. Energy error AE/E of two integrations of the Kepler prob- 
lem using the multiple timescale symplectic integrator with F c ^. The 
overall timestep tq = P/2000, Mj = 2, n = V2a, and rj/Vj+i = \/2. 



An alternative decomposition uses 
^p,i-l = < 



_ r / r 3 if r > rj, 

-/(^)^ 3 ifr m <r<n, (10) 
if r < rj+i, 



where /(x) = 2x 3 — 3x 2 + 1. Unlike the forces suggested by Skeel & Biesiadecki 
(1994), these forces have continuous first derivatives and decrease rapidly (or 
exactly) to zero at r <C r j (see Fig. 3). We shall not provide the details here, but 
we can understand why these properties are important from an analysis of H err . 
In Fig. 4 we show the energy error AE/E of two integrations using the multiple 
timescale symplectic integrator with F c ^. As expected, there is no secular drift 
in AE/E. Note also that, with the chosen integration parameters (r, oc rf), the 
maximum error is almost independent of the pericentric distance (which changes 
by 10 3 in the two cases shown). This again agrees with the expectation from an 
analysis of H err . 

One of the goals of this study is to develop a variable timestep integrator for 
solar system integrations. We have developed a second-order multiple timescale 
MVS integrator based on the algorithm described in this section (see Levison & 
Duncan 1994 for another approach). We are currently testing this integrator in 
detail. Initial results indicate that the integrator is fast and accurate and has 
all the desirable properties of the constant timestep symplectic integrators. 
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